% This file generates the impulse responses presented in Figures Figures 
% C2, C3, and C4 Appendix C of the paper. It is being called by 
% 'FIGs_C1_C2_C3_C4_Medium_Scale_Main.m'. This file runs the Dynare file 
% 'FIGs_C2_C3_C4_growth_base_model.mod'.

global y0 long drop rhoA replic iorder prune p_int peg sug dr ga alpha gi ...
       mcs rks xiw sigma pistar gammaw delta eta b beta ys chi phi ylsf klsf...
       eta beta b govs delta gup gi sigma wrsf wsf chi xiwf pistar gammaw prune

%% Options: 

fig_plot        =   0;
prune           =   1;          % set to 1 for pruning
growth_source   =   1;          % set to 0 for both A and I, 1 for just A, 2 for just I
third           =   0;          % set to 1 to do a third order approximation, 0 to do second

replic          =   50;         % number of replications when calculating irfs
drop            =   0;          % number of periods to drop when calculating irds
T               =   500;        % number of periods to simulate
burn            =   100;        % number of periods to burn-in after steady state
peg             =   8;          % number of periods to peg
long            =   22;         % length of IRFs


e11             =   clock;


alpha           =   1/3;
phi             =   0.001;
gad             =   1.0033;
gid             =   1.000;
thetag          =   0.5;        % parameter governing utility from government spending
gad             =   1;


iorder          =   1;

e1              =   clock;

%load estimated_params

psi1            =   1/((1-phi)*(1-alpha));
psi2            =   (alpha*(1-phi))/((1-phi)*(1-alpha));
check1          =   psi1*(phi-1 + alpha*(1-phi)) + 1;
check2          =   psi2*(phi-1 + alpha*(1-phi)) + alpha*(1-phi);

omegag          =   0.0;
beta            =   0.995;
theta           =   11;
sigma           =   11;
pistar          =   1.005;
pistar          =   1;
phi_fix         =   0.00;
% Desired output growth rate:
gupd            =   gad^(1/((1-phi_fix)*(1-alpha)))*gid^(alpha/(1-alpha)); 
gamma2          =   0.05;
delta           =   0.025;
xiw             =   0.66;
xip             =   0.66;
b               =   0.7;
chi             =   1;
kappa           =   4;
%rhoA            =   0.4;
gammap          =   0;
gammaw          =   0;
alphapi         =   1.5;
alphay          =   0.2;
rhoi            =   0.8;
rhoinvs         =   0.8;
rhob            =   0.6;
rhoL            =   0.8;
suA             =   0.01;
suL             =   0.01;
sug             =   0.01;
rhog            =   0.9;
sub             =   0.01;
sur             =   0.0025;
eta             =   6;
suinvs          =   0.04;
alphax          =   0;

if growth_source==0
    ga              =   gad^(1-phi);
    gi              =   gid;
    gup             =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
else
    if growth_source==1
        ga              =   gupd^((1-phi)*(1-alpha));
        gi              =   1;
        gup             =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
    else
        ga              =   1;
        gi              =   gupd^((1-alpha)/alpha);
        gup             =   ga^(1/((1-phi)*(1-alpha)))*gi^(alpha/(1-alpha));
    end
end
ga              =   1;
gxs             =   ga^(1/(1-alpha));

%% Compute steady states.

As              =   1;

ps              =   ((1-xip*pistar^((gammap-1)*(1-theta)))/(1-xip))^(1/(1-theta));
ss              =   (1-xip*pistar^(-theta*(gammap-1)))^(-1)*((1-xip)*ps^(-theta));
ints            =   pistar*gup/beta - 1;
mcs             =   ((theta-1)/theta)*((1-xip*pistar^((gammap-1)...
                    *(1-theta)))/(1-xip))^(1/(1-theta))*(1-xip*beta*pistar^((gammap-1)...
                    *(1-theta)))^(-1)*(1-xip*beta*pistar^(-theta*(gammap-1)));
rks             =   gup*gi/beta - (1-delta);
gamma1          =   rks;
ws              =   ((As*mcs*rks^(alpha*(phi-1)))/(phi^(-phi)*(1-phi)^(phi-1)...
                    *(alpha^(-alpha)*(1-alpha)^(alpha-1))^(1-phi)))^(1/((1-alpha)*(1-phi)));
wrs             =   ((1-xiw*gup^(sigma-1)*pistar^((gammaw-1)*(1-sigma)))/(1-xiw))^(1/(1-sigma))*ws;
kls             =   (alpha/(1-alpha))*gup*gi*(ws/rks);
gks             =   (phi/(1-phi))*(1/alpha)*gup^(-1)*gi^(-1)*rks;
xls             =   As*mcs*gks^(phi)*kls^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
xls2            =   As*(1/ss)*gks^(phi)*kls^(phi*(1-alpha) + alpha)*gup^(alpha*(phi-1))...
                    *gi^(alpha*(phi-1));
yls             =   xls - gks*kls;
cls             =   (1-omegag)*yls - (1-(1-delta)*gup^(-1)*gi^(-1))*kls;
ls              =   (((sigma-1)/sigma)*(cls)^(-1)*(1/eta)*(((gup-b) - beta*b...
                    *(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)))*(wrs/ws)^(sigma*chi)...
                    *wrs*(1-beta*xiw*gup^(sigma-1)*pistar^((sigma-1)...
                    *(1-gammaw)))^(-1)*(1-beta*xiw*gup^(sigma*(1+chi))...
                    *pistar^(sigma*(1+chi)*(1-gammaw))))^(1/(1+chi));
ys              =   yls*ls;
cs              =   cls*ls;
ks              =   kls*ls;
gs              =   gks*ks;
xs              =   xls*ls;
govs            =   omegag*ys;
invss           =   (1-(1-delta)*gup^(-1)*gi^(-1))*ks;
F               =   ((1-mcs*ss)/mcs)*xs;
% SS MUC:
lams            =   (1/cs)*((gup - b) - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)); 
x1s             =   (lams*xs*mcs)/(1-xip*beta*pistar^(-theta*(gammap-1)));
x2s             =   (lams*xs)/(1-xip*beta*pistar^((gammap-1)*(1-theta)));
f1s             =   (1-beta*xiw*gup^(sigma*(1+chi))*pistar^(sigma*(1+chi)...
                    *(1-gammaw)))^(-1)*(eta*(ws/wrs)^(sigma*(1+chi))*ls^(1+chi));
f2s             =   (1-beta*xiw*gup^(sigma-1)*pistar^((sigma-1)...
                    *(1-gammaw)))^(-1)*(lams*(ws/wrs)^(sigma)*ls);
% Check to make sure capital works out correctly:
k_check         =   ks - gi*gup*(alpha*(1-phi)*(mcs/rks)*(ss*xs + F)); 
% Check to make sure labor works out correctly:
l_check         =   ls - (1-alpha)*(1-phi)*(mcs/ws)*(ss*xs + F); 
% Check to make sure gamma works out:
gamma_check     =   gs - phi*mcs*(ss*xs + F); 
reset_price_check = ps - (theta/(theta-1))*x1s/x2s;
reset_wage_check =  wrs - (sigma/(sigma - 1))*f1s/f2s;
profit_check    =   xs - ws*ls - rks*ks*gup^(-1)*gi^(-1) - gs;
vws             =   (1-xiw*(gup^(-1)*pistar^(gammaw-1))^(-sigma*(1+chi)))^(-1)...
                    *((1-xiw)*(wrs/ws)^(-sigma*(1+chi)));

%% Find flexible price equilibrium.
xipf            =   0.00;
xiwf            =   0.00;

% Compute steady states. 
psf             =   ((1-xipf*pistar^((gammap-1)*(1-theta)))/(1-xipf))^(1/(1-theta));
ssf             =   (1-xipf*pistar^(-theta*(gammap-1)))^(-1)*((1-xipf)*psf^(-theta));
intsf           =   pistar*gup/beta - 1;
mcsf            =   ((theta-1)/theta)*((1-xipf*pistar^((gammap-1)...
                    *(1-theta)))/(1-xipf))^(1/(1-theta))*(1-xipf*beta...
                    *pistar^((gammap-1)*(1-theta)))^(-1)*(1-xipf*beta*pistar^(-theta*(gammap-1)));
rksf            =   gup*gi/beta - (1-delta);
wsf             =   ((mcsf*rksf^(alpha*(phi-1)))/(phi^(-phi)*(1-phi)^(phi-1)...
                    *(alpha^(-alpha)*(1-alpha)^(alpha-1))^(1-phi)))^(1/((1-alpha)*(1-phi)));
wrsf            =   ((1-xiwf*gup^(sigma-1)*pistar^((gammaw-1)*(1-sigma)))/(1-xiwf))^(1/(1-sigma))*wsf;
klsf            =   (alpha/(1-alpha))*gup*gi*(wsf/rksf);
gksf            =   (phi/(1-phi))*(1/alpha)*gup^(-1)*gi^(-1)*rksf;
xlsf            =   mcsf*gksf^(phi)*klsf^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
xls2f           =   (1/ss)*gksf^(phi)*klsf^(phi*(1-alpha) + alpha)...
                    *gup^(alpha*(phi-1))*gi^(alpha*(phi-1));
ylsf            =   xlsf - gksf*klsf;

% Numerically solve for C/L, G/L, and L.
options         =   optimset('TolX',1000e-25,'TolFun',1000e-25,...
                    'MaxFunEvals',1000e+25,'MaxIter',50,'LargeScale','off','display','iter');
guess           =   0.33;                   % Guess lsf.
[ll,fval]       =   fsolve(@FIGs_C1_C2_C3_C4_find_flex_price,guess,options);
lsf             =   ll;
clsf            =   ylsf - govs/lsf - (1-(1-delta)*gup^(-1)*gi^(-1))*klsf;
lsf_check       =   (((sigma-1)/sigma)*(clsf)^(-1)*(1/eta)*(((gup-b)...
                    - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))*(gup-b)))...
                    *(wrsf/wsf)^(sigma*chi)*wrsf*(1-beta*xiwf*gup^(sigma-1)...
                    *pistar^((sigma-1)*(1-gammaw)))^(-1)*(1-beta*xiwf...
                    *gup^(sigma*(1+chi))*pistar^(sigma*(1+chi)*(1-gammaw))))^(1/(1+chi));
ysf             =   ylsf*lsf;
csf             =   clsf*lsf;
ksf             =   klsf*lsf;
gsf             =   gksf*ksf;
xsf             =   xlsf*lsf;
govsf           =   omegag*ysf;
invssf          =   (1-(1-delta)*gup^(-1)*gi^(-1))*ksf;
Ff              =   ((1-mcsf*ssf)/mcsf)*xsf;
% SS MUC:
lamsf           =   (1/csf)*((gup - b) - beta*b*(1-b*gup^(-1)))/((1-b*gup^(-1))...
                    *(gup-b)); 
x1sf            =   (lamsf*xsf*mcsf)/(1-xipf*beta*pistar^(-theta*(gammap-1)));
x2sf            =   (lamsf*xsf)/(1-xipf*beta*pistar^((gammap-1)*(1-theta)));
f1sf            =   (1-beta*xiwf*gup^(sigma*(1+chi))*pistar^(sigma*(1+chi)...
                    *(1-gammaw)))^(-1)*(eta*(wsf/wrsf)^(sigma*(1+chi))*lsf^(1+chi));
f2sf            =   (1-beta*xiwf*gup^(sigma-1)*pistar^((sigma-1)...
                    *(1-gammaw)))^(-1)*(lamsf*(wsf/wrsf)^(sigma)*lsf);
% Check to make sure capital works out correctly:
k_checkf        =   ksf - gi*gup*(alpha*(1-phi)*(mcsf/rksf)*(ssf*xsf + Ff)); 
% Check to make sure labor works out correctly:
l_checkf        =   lsf - (1-alpha)*(1-phi)*(mcsf/wsf)*(ssf*xsf + Ff);
% Check to make sure gamma works out:
gamma_checkf    =   gsf - phi*mcsf*(ssf*xsf + Ff); 
reset_price_checkf = psf - (theta/(theta-1))*x1sf/x2sf;
reset_wage_checkf = wrsf - (sigma/(sigma - 1))*f1sf/f2sf;
profit_checkf   =   xsf - ws*lsf - rks*ksf*gup^(-1)*gi^(-1) - gsf;
vwsf            =   (1-xiwf*(gup^(-1)*pistar^(gammaw-1))^(-sigma*(1+chi)))^(-1)...
                     *((1-xiwf)*(wrsf/wsf)^(-sigma*(1+chi)));
gaps            =   ys/ysf;

%% Shocks, model solution, and variable positions. 


rhor            =   0.00;
css             =   cs;
sii             =   0.002;          % Standard deviation of anticipated shock.


save parameter_main_gov_peg beta b gamma1 gamma2 kappa delta pistar sigma...
     chi eta xiw gammaw alpha phi F theta xip gammap alphapi alphay rhoi...
     rhob rhoL rhoA rhor rhoinvs sub suL suA sur suinvs ints ys gup gi lams...
     ls ks ws wrs ps gs x1s x2s f1s f2s xs mcs ss As vws css lamsf csf wrsf...
     f1sf f2sf wsf lsf mcsf ssf xsf gsf psf x1sf x2sf ysf ksf invssf xipf xiwf...
     Ff gaps alphax omegag rhog sug thetag govs gxs

dynare FIGs_C2_C3_C4_growth_base_model noclearall nolog nostrict


dr              =   oo_.dr;
y0              =   oo_.dr.ys;

% Variable positions:
p_lam           =   1;
p_c             =   2;
p_rk            =   3;
p_z             =   4;
p_mu            =   5;
p_invs          =   6;
p_int           =   7;
p_infl          =   8;
p_wr            =   9;
p_f1            =   10;
p_f2            =   11;
p_khat          =   12;
p_v             =   13;
p_s             =   14;
p_L             =   15;
p_pr            =   16;
p_x1            =   17;
p_x2            =   18;
p_w             =   19;
p_y             =   20;
p_k             =   21;
p_eb            =   22;
p_eL            =   23;
p_A             =   24;
p_er            =   25;
p_einvs         =   26;
p_dy            =   27;
p_dc            =   28;
p_dI            =   29;
p_dL            =   30;
p_dw            =   31;

guess           =   zeros(1,peg);
[xx,fval]       =   fsolve(@FIGs_C1_C2_C3_C4_find_peg_shocks,guess);
e1              =   [0 0 0.01 xx(1) 0 0 xx(2) xx(3) xx(4) xx(5) xx(6) xx(7) xx(8)]';

%% Impulse Response Functions. 

% Modification of Dynare's IRF command:
irf_ss          =   FIGs_C1_C2_C3_C4_irf_dynare(dr, e1, long+1, drop, replic, iorder); 
r_irf_ss        =   zeros(1,long);

for jx = 1:long
    r_irf_ss(jx)    =   irf_ss(p_int,jx) - irf_ss(p_infl,jx+1);
end

% Now do it without the ZLB:
xx              =   zeros(1,peg);
e2              =   [0 0 0.01 xx(1) 0 0 xx(2) xx(3) xx(4) xx(5) xx(6) xx(7) xx(8)]';
irf2            =   FIGs_C1_C2_C3_C4_irf_dynare(dr, e2, long+1, drop, replic, iorder); 
r_irf2          =   zeros(1,long);

for jx = 1:long
    r_irf2(jx)      =   irf2(p_int,jx) - irf2(p_infl,jx+1);
end

e22             =   clock;

elapsed         =   etime(e22,e11)/60